Load Data
dataset <- read.delim("raw_data/FigureS5J.txt", stringsAsFactors = FALSE)
dataset$genotype <- gsub(" ", "", dataset$genotype )
dataset$genotype <- factor(dataset$genotype)
dataset$Experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each=length(unique(dataset$genotype))))
dataset$UID <- factor(paste(dataset$Experiment, dataset$genotype))
# wide format
kable(dataset, row.names = F)
| WT |
1160 |
810 |
567 |
151 |
exp1 |
exp1 WT |
| YFP-ALC1 |
1005 |
792 |
585 |
192 |
exp1 |
exp1 YFP-ALC1 |
| WT |
1056 |
838 |
378 |
258 |
exp2 |
exp2 WT |
| YFP-ALC1 |
2132 |
1742 |
983 |
399 |
exp2 |
exp2 YFP-ALC1 |
| WT |
1604 |
1036 |
558 |
184 |
exp3 |
exp3 WT |
| YFP-ALC1 |
1714 |
1325 |
985 |
311 |
exp3 |
exp3 YFP-ALC1 |
library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "Treatment", value.name = "Counts")
dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT")
dataset$Bleomycin <- gsub("NT","1",dataset$Treatment)
dataset$Bleomycin <- gsub("bleomycin_|uM","",dataset$Bleomycin)
dataset$Bleomycin <- log10(as.integer(dataset$Bleomycin))
dataset$Offset <- NA
for(uid in levels(dataset$UID)){
dataset$Offset[dataset$UID == uid] <- mean(dataset$Counts[dataset$UID == uid])
}
dataset$NormCounts <- dataset$Counts / dataset$Offset
dataset$Offset2 <- NA
for(gid in levels(dataset$genotype)){
dataset$Offset2[dataset$genotype == gid] <- mean(dataset$NormCounts[dataset$genotype == gid & dataset$Bleomycin == 0])
}
dataset$NormCounts2 <- dataset$NormCounts / dataset$Offset2
# long format
kable(dataset, row.names = F)
| WT |
exp1 |
exp1 WT |
NT |
1160 |
0.00000 |
672.00 |
1.7261905 |
1.764286 |
0.9784074 |
| YFP-ALC1 |
exp1 |
exp1 YFP-ALC1 |
NT |
1005 |
0.00000 |
643.50 |
1.5617716 |
1.588615 |
0.9831029 |
| WT |
exp2 |
exp2 WT |
NT |
1056 |
0.00000 |
632.50 |
1.6695652 |
1.764286 |
0.9463121 |
| YFP-ALC1 |
exp2 |
exp2 YFP-ALC1 |
NT |
2132 |
0.00000 |
1314.00 |
1.6225266 |
1.588615 |
1.0213469 |
| WT |
exp3 |
exp3 WT |
NT |
1604 |
0.00000 |
845.50 |
1.8971023 |
1.764286 |
1.0752805 |
| YFP-ALC1 |
exp3 |
exp3 YFP-ALC1 |
NT |
1714 |
0.00000 |
1083.75 |
1.5815456 |
1.588615 |
0.9955502 |
| WT |
exp1 |
exp1 WT |
bleomycin_5uM |
810 |
0.69897 |
672.00 |
1.2053571 |
1.764286 |
0.6831983 |
| YFP-ALC1 |
exp1 |
exp1 YFP-ALC1 |
bleomycin_5uM |
792 |
0.69897 |
643.50 |
1.2307692 |
1.588615 |
0.7747438 |
| WT |
exp2 |
exp2 WT |
bleomycin_5uM |
838 |
0.69897 |
632.50 |
1.3249012 |
1.764286 |
0.7509560 |
| YFP-ALC1 |
exp2 |
exp2 YFP-ALC1 |
bleomycin_5uM |
1742 |
0.69897 |
1314.00 |
1.3257230 |
1.588615 |
0.8345152 |
| WT |
exp3 |
exp3 WT |
bleomycin_5uM |
1036 |
0.69897 |
845.50 |
1.2253105 |
1.764286 |
0.6945078 |
| YFP-ALC1 |
exp3 |
exp3 YFP-ALC1 |
bleomycin_5uM |
1325 |
0.69897 |
1083.75 |
1.2226067 |
1.588615 |
0.7696056 |
| WT |
exp1 |
exp1 WT |
bleomycin_10uM |
567 |
1.00000 |
672.00 |
0.8437500 |
1.764286 |
0.4782388 |
| YFP-ALC1 |
exp1 |
exp1 YFP-ALC1 |
bleomycin_10uM |
585 |
1.00000 |
643.50 |
0.9090909 |
1.588615 |
0.5722539 |
| WT |
exp2 |
exp2 WT |
bleomycin_10uM |
378 |
1.00000 |
632.50 |
0.5976285 |
1.764286 |
0.3387367 |
| YFP-ALC1 |
exp2 |
exp2 YFP-ALC1 |
bleomycin_10uM |
983 |
1.00000 |
1314.00 |
0.7480974 |
1.588615 |
0.4709118 |
| WT |
exp3 |
exp3 WT |
bleomycin_10uM |
558 |
1.00000 |
845.50 |
0.6599645 |
1.764286 |
0.3740689 |
| YFP-ALC1 |
exp3 |
exp3 YFP-ALC1 |
bleomycin_10uM |
985 |
1.00000 |
1083.75 |
0.9088812 |
1.588615 |
0.5721219 |
| WT |
exp1 |
exp1 WT |
bleomycin_50uM |
151 |
1.69897 |
672.00 |
0.2247024 |
1.764286 |
0.1273617 |
| YFP-ALC1 |
exp1 |
exp1 YFP-ALC1 |
bleomycin_50uM |
192 |
1.69897 |
643.50 |
0.2983683 |
1.588615 |
0.1878167 |
| WT |
exp2 |
exp2 WT |
bleomycin_50uM |
258 |
1.69897 |
632.50 |
0.4079051 |
1.764286 |
0.2312013 |
| YFP-ALC1 |
exp2 |
exp2 YFP-ALC1 |
bleomycin_50uM |
399 |
1.69897 |
1314.00 |
0.3036530 |
1.588615 |
0.1911433 |
| WT |
exp3 |
exp3 WT |
bleomycin_50uM |
184 |
1.69897 |
845.50 |
0.2176227 |
1.764286 |
0.1233489 |
| YFP-ALC1 |
exp3 |
exp3 YFP-ALC1 |
bleomycin_50uM |
311 |
1.69897 |
1083.75 |
0.2869666 |
1.588615 |
0.1806395 |
Plot Data
library(ggplot2)
# raw data
ggplot(dataset, aes(x=Bleomycin, y=Counts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=genotype)) +
geom_point(aes(colour=genotype, shape=Experiment), size=2) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
scale_shape_manual(values=15:20) +
scale_color_manual(values=c("#000000","#808000"))

# NormCounts Linear
ggplot(dataset, aes(x=Bleomycin, y=NormCounts, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=Bleomycin, y=NormCounts2, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
scale_color_manual(values=c("#000000","#808000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=Bleomycin, y=NormCounts, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)")+
scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=Bleomycin, y=NormCounts2, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
scale_color_manual(values=c("#000000","#808000"))

# NormCounts Cubic
ggplot(dataset, aes(x=Bleomycin, y=NormCounts, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)")+
scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=Bleomycin, y=NormCounts2, color=genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=genotype), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
scale_color_manual(values=c("#000000","#808000"))

library(Cairo)
cairo_pdf("FigureS5J.pdf", width = 5, height = 4, family = "Arial")
ggplot(dataset, aes(x=Bleomycin, y=NormCounts2)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank()) +
geom_point(aes(colour = genotype)) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=TRUE, aes(colour = genotype), fill='#CCCCCC') +
#facet_grid(. ~ genotype) +
xlab(label = "Bleomycin (log10 µM)") +
ylab(label = "Normalized Counts") +
scale_color_manual(values=c("#000000","#808000"))
dev.off()
## quartz_off_screen
## 2
Models
library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)
Linear formula
fit1 <- lm(Counts ~ Experiment + Bleomycin*genotype, data = dataset)
print(summary(fit1))
##
## Call:
## lm(formula = Counts ~ Experiment + Bleomycin * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -476.82 -107.86 -15.27 130.19 500.52
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1064.5 141.8 7.509 5.97e-07 ***
## Experimentexp2 315.5 123.7 2.550 0.0201 *
## Experimentexp3 306.9 123.7 2.481 0.0232 *
## Bleomycin -653.7 117.1 -5.583 2.68e-05 ***
## genotypeYFP-ALC1 417.3 173.2 2.410 0.0269 *
## Bleomycin:genotypeYFP-ALC1 -141.6 165.6 -0.855 0.4038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 247.4 on 18 degrees of freedom
## Multiple R-squared: 0.8398, Adjusted R-squared: 0.7953
## F-statistic: 18.88 on 5 and 18 DF, p-value: 1.355e-06
cat("AIC: ", AIC(fit1))
## AIC: 339.7383
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ Bleomycin*genotype, data = dataset)
print(summary(fit2))
##
## Call:
## lm(formula = NormCounts ~ Bleomycin * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26678 -0.05879 -0.01889 0.09365 0.20857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.76526 0.06366 27.729 < 2e-16 ***
## Bleomycin -0.90086 0.06087 -14.799 3.08e-12 ***
## genotypeYFP-ALC1 -0.10406 0.09003 -1.156 0.261
## Bleomycin:genotypeYFP-ALC1 0.12249 0.08608 1.423 0.170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1286 on 20 degrees of freedom
## Multiple R-squared: 0.9503, Adjusted R-squared: 0.9429
## F-statistic: 127.5 on 3 and 20 DF, p-value: 3.315e-13
cat("AIC: ", AIC(fit2))
## AIC: -24.70436
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ Bleomycin*genotype, data = dataset)
print(summary(fit3))
##
## Call:
## lm(formula = NormCounts2 ~ Bleomycin * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15121 -0.03700 -0.01070 0.05473 0.13129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.00055 0.03746 26.709 < 2e-16 ***
## Bleomycin -0.51061 0.03582 -14.255 6.13e-12 ***
## genotypeYFP-ALC1 0.04514 0.05298 0.852 0.404
## Bleomycin:genotypeYFP-ALC1 0.02064 0.05066 0.408 0.688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07569 on 20 degrees of freedom
## Multiple R-squared: 0.9517, Adjusted R-squared: 0.9445
## F-statistic: 131.5 on 3 and 20 DF, p-value: 2.479e-13
cat("AIC: ", AIC(fit3))
## AIC: -50.15771
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ Bleomycin*genotype + (1|UID), data = dataset)
print(summary(fit4))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ Bleomycin * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 282.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.91400 -0.53151 -0.01178 0.60395 1.83213
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 55274 235.1
## Residual 36717 191.6
## Number of obs: 24, groups: UID, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1271.943 165.583 6.393 7.682 0.000187 ***
## Bleomycin -653.663 90.676 16.000 -7.209 2.08e-06 ***
## genotypeYFP-ALC1 417.333 234.169 6.393 1.782 0.121961
## Bleomycin:genotypeYFP-ALC1 -141.556 128.235 16.000 -1.104 0.285973
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Blmycn gYFP-A
## Bleomycin -0.465
## gntYFP-ALC1 -0.707 0.329
## Bl:YFP-ALC1 0.329 -0.707 -0.465
cat("AIC: ", AIC(fit4))
## AIC: 294.7363
simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula
fit5 <- lm(Counts ~ Experiment + poly(Bleomycin,2)*genotype, data = dataset)
print(summary(fit5))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(Bleomycin, 2) * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -421.90 -106.95 -25.02 149.03 445.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 509.21 105.15 4.842 0.00018 ***
## Experimentexp2 315.50 128.79 2.450 0.02618 *
## Experimentexp3 306.88 128.79 2.383 0.02993 *
## poly(Bleomycin, 2)1 -1953.50 364.27 -5.363 6.34e-05 ***
## poly(Bleomycin, 2)2 92.26 364.27 0.253 0.80327
## genotypeYFP-ALC1 297.08 105.15 2.825 0.01219 *
## poly(Bleomycin, 2)1:genotypeYFP-ALC1 -423.04 515.15 -0.821 0.42360
## poly(Bleomycin, 2)2:genotypeYFP-ALC1 -361.30 515.15 -0.701 0.49316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 257.6 on 16 degrees of freedom
## Multiple R-squared: 0.8457, Adjusted R-squared: 0.7782
## F-statistic: 12.53 on 7 and 16 DF, p-value: 1.973e-05
cat("AIC: ", AIC(fit5))
## AIC: 342.8408
simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(Bleomycin,2)*genotype, data = dataset)
print(summary(fit6))
##
## Call:
## lm(formula = NormCounts ~ poly(Bleomycin, 2) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.242931 -0.041198 0.004446 0.065468 0.213157
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+00 3.629e-02 27.559 3.58e-16 ***
## poly(Bleomycin, 2)1 -2.692e+00 1.778e-01 -15.145 1.10e-11 ***
## poly(Bleomycin, 2)2 1.168e-01 1.778e-01 0.657 0.519
## genotypeYFP-ALC1 4.392e-17 5.132e-02 0.000 1.000
## poly(Bleomycin, 2)1:genotypeYFP-ALC1 3.661e-01 2.514e-01 1.456 0.163
## poly(Bleomycin, 2)2:genotypeYFP-ALC1 -3.987e-01 2.514e-01 -1.586 0.130
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1257 on 18 degrees of freedom
## Multiple R-squared: 0.9573, Adjusted R-squared: 0.9454
## F-statistic: 80.71 on 5 and 18 DF, p-value: 1.097e-11
cat("AIC: ", AIC(fit6))
## AIC: -24.34165
simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(Bleomycin,2)*genotype, data = dataset)
print(summary(fit7))
##
## Call:
## lm(formula = NormCounts2 ~ poly(Bleomycin, 2) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.137694 -0.024008 0.002698 0.039744 0.120818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.56680 0.02115 26.794 5.88e-16 ***
## poly(Bleomycin, 2)1 -1.52597 0.10363 -14.724 1.76e-11 ***
## poly(Bleomycin, 2)2 0.06622 0.10363 0.639 0.5309
## genotypeYFP-ALC1 0.06268 0.02992 2.095 0.0506 .
## poly(Bleomycin, 2)1:genotypeYFP-ALC1 0.06170 0.14656 0.421 0.6788
## poly(Bleomycin, 2)2:genotypeYFP-ALC1 -0.24363 0.14656 -1.662 0.1138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07328 on 18 degrees of freedom
## Multiple R-squared: 0.9593, Adjusted R-squared: 0.948
## F-statistic: 84.83 on 5 and 18 DF, p-value: 7.161e-12
cat("AIC: ", AIC(fit7))
## AIC: -50.24129
simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(Bleomycin,2)*genotype + (1|UID), data = dataset)
print(summary(fit8))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(Bleomycin, 2) * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 251.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.59469 -0.57627 -0.04532 0.59448 1.51209
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 54685 233.8
## Residual 39074 197.7
## Number of obs: 24, groups: UID, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 716.67 146.57 4.00 4.889
## poly(Bleomycin, 2)1 -1953.49 279.55 14.00 -6.988
## poly(Bleomycin, 2)2 92.26 279.55 14.00 0.330
## genotypeYFP-ALC1 297.08 207.29 4.00 1.433
## poly(Bleomycin, 2)1:genotypeYFP-ALC1 -423.05 395.34 14.00 -1.070
## poly(Bleomycin, 2)2:genotypeYFP-ALC1 -361.30 395.34 14.00 -0.914
## Pr(>|t|)
## (Intercept) 0.00811 **
## poly(Bleomycin, 2)1 6.37e-06 ***
## poly(Bleomycin, 2)2 0.74625
## genotypeYFP-ALC1 0.22509
## poly(Bleomycin, 2)1:genotypeYFP-ALC1 0.30269
## poly(Bleomycin, 2)2:genotypeYFP-ALC1 0.37624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(B,2)1 pl(B,2)2 gYFP-A p(B,2)1:
## ply(Blm,2)1 0.000
## ply(Blm,2)2 0.000 0.000
## gntYFP-ALC1 -0.707 0.000 0.000
## p(B,2)1:YFP 0.000 -0.707 0.000 0.000
## p(B,2)2:YFP 0.000 0.000 -0.707 0.000 0.000
cat("AIC: ", AIC(fit8))
## AIC: 267.1439
simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula
fit9 <- lm(Counts ~ Experiment + poly(Bleomycin,3)*genotype, data = dataset)
print(summary(fit9))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(Bleomycin, 3) * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -404.54 -95.08 -6.06 104.79 406.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 509.208 105.905 4.808 0.000278 ***
## Experimentexp2 315.500 129.707 2.432 0.029009 *
## Experimentexp3 306.875 129.707 2.366 0.032948 *
## poly(Bleomycin, 3)1 -1953.495 366.867 -5.325 0.000107 ***
## poly(Bleomycin, 3)2 92.264 366.867 0.251 0.805088
## poly(Bleomycin, 3)3 346.343 366.867 0.944 0.361149
## genotypeYFP-ALC1 297.083 105.905 2.805 0.014036 *
## poly(Bleomycin, 3)1:genotypeYFP-ALC1 -423.045 518.828 -0.815 0.428506
## poly(Bleomycin, 3)2:genotypeYFP-ALC1 -361.300 518.828 -0.696 0.497592
## poly(Bleomycin, 3)3:genotypeYFP-ALC1 -1.664 518.828 -0.003 0.997487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 259.4 on 14 degrees of freedom
## Multiple R-squared: 0.8631, Adjusted R-squared: 0.775
## F-statistic: 9.804 on 9 and 14 DF, p-value: 0.0001174
cat("AIC: ", AIC(fit9))
## AIC: 343.9775
simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(Bleomycin,3)*genotype, data = dataset)
print(summary(fit10))
##
## Call:
## lm(formula = NormCounts ~ poly(Bleomycin, 3) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10726 -0.04199 -0.01795 0.05358 0.14330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+00 2.481e-02 40.307 < 2e-16 ***
## poly(Bleomycin, 3)1 -2.692e+00 1.215e-01 -22.151 1.97e-13 ***
## poly(Bleomycin, 3)2 1.168e-01 1.215e-01 0.961 0.350736
## poly(Bleomycin, 3)3 4.929e-01 1.215e-01 4.056 0.000918 ***
## genotypeYFP-ALC1 8.710e-17 3.509e-02 0.000 1.000000
## poly(Bleomycin, 3)1:genotypeYFP-ALC1 3.661e-01 1.719e-01 2.130 0.049063 *
## poly(Bleomycin, 3)2:genotypeYFP-ALC1 -3.987e-01 1.719e-01 -2.319 0.033931 *
## poly(Bleomycin, 3)3:genotypeYFP-ALC1 -1.938e-01 1.719e-01 -1.128 0.276093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08594 on 16 degrees of freedom
## Multiple R-squared: 0.9823, Adjusted R-squared: 0.9745
## F-statistic: 126.5 on 7 and 16 DF, p-value: 8.384e-13
cat("AIC: ", AIC(fit10))
## AIC: -41.41727
simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(Bleomycin,3)*genotype, data = dataset)
print(summary(fit11))
##
## Call:
## lm(formula = NormCounts2 ~ poly(Bleomycin, 3) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06752 -0.02410 -0.01047 0.03373 0.08122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.56680 0.01442 39.314 < 2e-16 ***
## poly(Bleomycin, 3)1 -1.52597 0.07063 -21.605 2.9e-13 ***
## poly(Bleomycin, 3)2 0.06622 0.07063 0.938 0.36241
## poly(Bleomycin, 3)3 0.27939 0.07063 3.956 0.00113 **
## genotypeYFP-ALC1 0.06268 0.02039 3.074 0.00726 **
## poly(Bleomycin, 3)1:genotypeYFP-ALC1 0.06170 0.09989 0.618 0.54549
## poly(Bleomycin, 3)2:genotypeYFP-ALC1 -0.24363 0.09989 -2.439 0.02675 *
## poly(Bleomycin, 3)3:genotypeYFP-ALC1 -0.09112 0.09989 -0.912 0.37520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04994 on 16 degrees of freedom
## Multiple R-squared: 0.9832, Adjusted R-squared: 0.9758
## F-statistic: 133.7 on 7 and 16 DF, p-value: 5.446e-13
cat("AIC: ", AIC(fit11))
## AIC: -67.4721
simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(Bleomycin,3)*genotype + (1|UID), data = dataset)
print(summary(fit12))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(Bleomycin, 3) * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 221.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55170 -0.46246 -0.07901 0.43109 1.35742
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 55544 235.7
## Residual 35638 188.8
## Number of obs: 24, groups: UID, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 716.667 146.575 4.000 4.889
## poly(Bleomycin, 3)1 -1953.495 266.975 12.000 -7.317
## poly(Bleomycin, 3)2 92.264 266.975 12.000 0.346
## poly(Bleomycin, 3)3 346.343 266.975 12.000 1.297
## genotypeYFP-ALC1 297.083 207.289 4.000 1.433
## poly(Bleomycin, 3)1:genotypeYFP-ALC1 -423.045 377.560 12.000 -1.120
## poly(Bleomycin, 3)2:genotypeYFP-ALC1 -361.300 377.560 12.000 -0.957
## poly(Bleomycin, 3)3:genotypeYFP-ALC1 -1.664 377.560 12.000 -0.004
## Pr(>|t|)
## (Intercept) 0.00811 **
## poly(Bleomycin, 3)1 9.26e-06 ***
## poly(Bleomycin, 3)2 0.73563
## poly(Bleomycin, 3)3 0.21892
## genotypeYFP-ALC1 0.22509
## poly(Bleomycin, 3)1:genotypeYFP-ALC1 0.28444
## poly(Bleomycin, 3)2:genotypeYFP-ALC1 0.35748
## poly(Bleomycin, 3)3:genotypeYFP-ALC1 0.99656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(B,3)1 pl(B,3)2 pl(B,3)3 gYFP-A p(B,3)1: p(B,3)2:
## ply(Blm,3)1 0.000
## ply(Blm,3)2 0.000 0.000
## ply(Blm,3)3 0.000 0.000 0.000
## gntYFP-ALC1 -0.707 0.000 0.000 0.000
## p(B,3)1:YFP 0.000 -0.707 0.000 0.000 0.000
## p(B,3)2:YFP 0.000 0.000 -0.707 0.000 0.000 0.000
## p(B,3)3:YFP 0.000 0.000 0.000 -0.707 0.000 0.000 0.000
cat("AIC: ", AIC(fit12))
## AIC: 241.8309
simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Compare Results
ICtab(fit1,fit2,fit3,fit4,
fit5,fit6,fit7,fit8,
fit9,fit10,fit11,fit12,
base=T)
## AIC dAIC df
## fit11 -67.5 0.0 9
## fit7 -50.2 17.2 7
## fit3 -50.2 17.3 5
## fit10 -41.4 26.1 9
## fit2 -24.7 42.8 5
## fit6 -24.3 43.1 7
## fit12 241.8 309.3 10
## fit8 267.1 334.6 8
## fit4 294.7 362.2 6
## fit1 339.7 407.2 7
## fit5 342.8 410.3 9
## fit9 344.0 411.4 11
Final Result
fit <- fit11
output <- coef(summary(fit))
output <- output[grep("Bleomycin", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
# suggested result table
kable(output, row.names = T)
| Bleomycin1 in WT |
-1.5259664 |
0.0706299 |
-21.6050932 |
0.0000000 |
| Bleomycin2 in WT |
0.0662200 |
0.0706299 |
0.9375630 |
0.3624070 |
| Bleomycin3 in WT |
0.2793890 |
0.0706299 |
3.9556733 |
0.0011329 |
| Bleomycin1: WT vs. YFP-ALC1 |
0.0616959 |
0.0998858 |
0.6176642 |
0.5454882 |
| Bleomycin2: WT vs. YFP-ALC1 |
-0.2436254 |
0.0998858 |
-2.4390393 |
0.0267529 |
| Bleomycin3: WT vs. YFP-ALC1 |
-0.0911161 |
0.0998858 |
-0.9122025 |
0.3752023 |
write.table(output, file = "FigureS5J_Stats.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
Anova
fit11a <- lm(NormCounts2 ~ poly(Bleomycin,3)*genotype, data = dataset)
fit11b <- lm(NormCounts2 ~ poly(Bleomycin,3)+genotype, data = dataset)
# anova table
anova(fit11a, fit11b)
## Analysis of Variance Table
##
## Model 1: NormCounts2 ~ poly(Bleomycin, 3) * genotype
## Model 2: NormCounts2 ~ poly(Bleomycin, 3) + genotype
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 16 0.039909
## 2 19 0.057774 -3 -0.017865 2.3875 0.1071